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Foreword 


Efficient  and  robust  ocean  models  are  an  essential  component 
of  any  ocean  prediction  system.  ITie  ocean  model  described  here 
is  part  of  several  prediction  products  currently  under  evaluation  by 
the  Navy.  In  conjunction  with  this  report,  the  source  code  for  the 
ocean  model  is  being  made  available  to  the  ocean  modeling  com¬ 
munity.  The  Naval  Oceanographic  and  Atmospheric  Research 
Laboratory  hopes  that  it  will  be  joined  by  other  ocean  models  to 
form  a  community  resource  library,  as  recommended  by  the  recent 
Workshop  on  Ocean  Prediction. 

W.  B.  Moseley  L.  R.  Elliott,  Commander,  USN 

Technical  Director  Commanding  Officer 


Executive  Summary  This  report  is  a  users  guide  to  the  Navy’s  hydrodynamic  (isopycnal) 

nonlinear,  primitive  equation,  layered  ocean  circulation  model.  The  model 
retains  the  free  surface  and  uses  a  semi-implicit  time  scheme  that  treats 
all  gravity  waves  implicitly.  It  can  handle  full-scale  bottom  topography, 
provided  it  is  confined  to  the  lowest  layer,  and  an  arbitrary  coastline 
geometry. 

The  model  has  been  in  use  at  the  Naval  Oceanographic  and  Atmo¬ 
spheric  Research  Laboratory  for  more  than  10  years  for  simulations  of 
the  ocean  circulation  in  the  Gulf  of  Mexico,  the  Caribbean  Sea,  the 
Alboran  Sea,  the  western  Mediterranean  Sea,  and  the  global  oceans.  In 
conjunction  with  the  issuance  of  this  report,  the  model  code  is  being 
made  available  to  the  ocean  modeling  community. 

The  vertically  integrated  equations  of  motion  and  their  finite  difference 
discretization  on  a  C-grid  is  presented,  as  is  a  description  on  the  semi- 
implicit  time  scheme,  the  boundary  conditions,  and  the  external  forcing. 

The  model  code  contains  internal  documentation  that  fully  describes 
the  user-specified  model  parameters  and  data  sets.  This  report  also 
contains  general  information  about  how  to  use  the  model,  in  particular, 
how  to  set  it  up  for  a  new  ocean  region  and  how  to  port  it  to  a  new 
computer  system. 


This  work  was  funded  by  the  Naval  Air  Development  Center,  Program  Acknowledgments 
Element  0602435N.  George  Hebum,  Harley  Hurlburt,  Dana  Thompson, 
and  many  other  past  and  present  members  of  NOARL’s  Ocean 
Dynamics  and  Prediction  Branch  contributed  to  the  development  of  the 
ocean  model  described  in  this  report. 

The  mention  of  commercial  products  or  the  use  of  company  names 
does  not  in  any  way  imply  endorsement  by  the  U.S.  Navy  or  NOARL. 
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1.0  Introduction  The  layered  ocean  circulation  model  has  been  under  development 

since  1976  by  the  Naval  Oceanographic  and  Atmospheric  Research 
Laboratory  (NOARL).  TTie  original  version,  written  by  Hurlburt  and 
Thompson  (1980),  was  a  one-  or  two-layer  finite  depth,  or  one- 
active-layer  reduced-gravity,  hydrodynamic  primitive  equation  ocean 
model  on  a  beta-plane.  This  version  retained  the  free  surface  and  used 
a  semi-implicit  time  scheme,  was  formulated  on  the  C-grid,  and  could 
handle  only  rectangular  or  simple  nonrectangular  regions.  The  Helmholtz 
equations  from  the  semi-implicit  scheme  were  solved  using  a  version  of 
the  Hockney  (1965)  method,  as  implemented  by  Moore  (1989).  The 
capacitance  matrix  technique  (Buzbee  et  al.,  1971)  was  used  for 
nonrectangular  regions. 

The  model  was  first  used  at  NOARL  in  a  major  study  of  the  Gulf  of 
Mexico  (Hurlburt  and  Thompson,  1980;  1982),  and  a  description  of  the 
model  as  it  then  existed  was  included  as  an  appendix  to  the  1980  paper. 
Since  that  time  the  capability  of  the  model  has  been  greatly  expanded 
and  the  computer  code  completely  rewritten.  Hurlburt  and  Thompson 
(1980)  remain  the  appropriate  refereed  reference  on  the  basic  model 
design.  Their  article  includes  many  details  that  will  not  be  repeated 
here  and  should  be  used  in  addition  to  this  report. 

The  major  features  of  the  Navy  layered  ocean  circulation  model  are 
listed. 

•  Layers  in  the  vertical. 

•  Nonlinear  primitive  equations. 

•  Free  surface. 

•  Semi-implicit  time  scheme. 

•  Hydrodynamic,  i.e.,  isopycnal. 

•  Full-scale  bottom  topography  in  lowest  layer. 

•  Arbitrary  coastline  geometry 

Upgrades  to  the  original  model  include  the  following. 

•  The  model  region  can  be  any  closed  geometry;  i.e.,  any  coastline 
resolvable  by  the  grid  spacing  can  be  used.  There  is  no  limit  on  the 
number  of  islands  in  the  region. 

•  The  model  grid  is  uniform,  but  it  can  be  on  a  beta-plane,  on 
an  f-plane,  or  on  the  surface  of  an  earth-sized  sphere.  The  spherical 
geometry  version  was  formulated  by  Hebum  (NOARL,  personal 
communication). 

•  The  model  can  have  any  number  of  layers. 

•  The  external  mode  Helmholtz  equation  from  the  semi-implicit  lime 
scheme  is  solved  by  a  standard  FORTRAN  77  package  that  u.scs  a  fast 
direct  solver  for  rectangular  regions  (Swarztrauber,  1977)  and  the 
capacitance  matrix  technique  for  nonrectangular  regions.  Moore  (1989) 
wrote  the  fast  Fourier  transform  component  of  the  fast  direct  solver. 
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•  The  internal  mode  Helmholtz  equations  are  solved  by  a  standard 
FORTRAN  77  implementation  of  the  red-black  successive  over-relaxation 
iterative  scheme  (Hageman  and  Young,  1981). 

•  Versions  of  the  reduced  gravity  models  that  use  an  explicit  time 
step  are  available. 

•  Open  outflow  boundaries  now  use  a  variant  of  the  Orlanski  (1976) 
radiation  condition,  developed  by  Hebum  (1986). 

•  Any  number  of  open-boundary  sections  can  be  placed  along  any 
zonal  or  meridional  sections  of  “coastline,”  and  each  layer  can  have  its 
own  set  of  open  boundary  sections. 

•  “Hydromixing”  may  be  used  to  control  the  surfacing  of  layer  inter¬ 
faces,  using  a  formulation  by  Thompson  and  Wallcraft  (personal 
communication).  It  is  identical  to  mixing  in  a  thermodynamic  layered 
model,  except  that  layer  density  is  not  changed. 

The  model  has  been  optimized  to  study  and  predict  mesoscale  oceanic 
features  over  relatively  large  regions.  It  has  been  used  at  NOARL  to 
study  the  ocean  circulation  of  the  following  regions. 

•  The  Gulf  of  Mexico  at  grid  resolutions  of  0.2°  and  0.1°  (Hurlburt 
and  Thompson,  1980  and  1982;  Hurlburt,  1984  and  1986;  Hurlburt 
et  al.,  1990;  Kindle,  1986;  Thompson,  1986;  Wallcraft,  1985  and  1986). 

•  The  Caribbean  (Hebum  et  al.,  1982;  Kinder  et  al.,  1985). 

•  The  AJboran  Sea  (Preller  and  Hurlburt,  1982;  Parrilla  et  al.,  1986; 
Preller,  1986). 

•  The  Western  Mediterranean  at  7.5  by  5  km  (Hebum,  1986  and 
1988;  Hebum  and  La  Violette,  1990). 

•  The  Gulf  Stream  region  at  0.2°  and  0.125°  resolution  (Hurlburt  and 
Thompson,  1984;  Thompson  and  Schmitz,  1989). 

•  The  Indian  Ocean  at  0.4°  (Kindle  and  Thompson,  19o9). 

•The  Pacific  north  of  20°S  at  0.25°  and  0.125°  (Hurlburt  et  al., 
1991). 

•  Tlie  global  ocean  at  0.5°  resolution  (Hurlburt  et  al.,  1989;  Kindle 
et  al.,  1989;  Murray  ct  al.,  1989  and  1990). 

Layers  are  an  efficient  way  of  modeling  the  ocean’s  vertical  structure. 
Successful  simulations  are  possible  with  two  or  three  appropriately  located 
Lagrangian  layers,  while  at  least  20  vertical  levels  are  typically  used  in 
level  model  simulations  (Semtner  and  Chervin,  1988).  Moreover,  the 
model’s  semi-implicit  time  scheme  treats  all  gravity  waves  implicitly, 
which  allows  a  longer  time  step  than  primitive  equation  models  that  use 
a  rigid  lid  or  time  splitting.  The  difference  in  time  step  depends  on  the 
region  being  simulated,  but  typically  the  semi-implicit  time  step  is  three 
times  longer  than  that  required  by  schemes  that  treat  internal  gravity 
waves  explicitly.  The  model  is  therefore  at  least  100  times  more  efficient, 
in  terms  of  central  processing  time,  than  existing  level-type  primitive 
equation  ocean  models  with  20  or  more  levels,  and  it  also  requires 
significantly  less  computer  memory.  The  major  disadvantage  of  this 
layer  formulation  is  that  bottom  topography  is  confined  to  the 
lowest  layer,  which  makes  the  model  less  u,seful  in  coastal  regions. 
However,  it  has  proved  quite  successful  at  modeling  seamounts,  e.g., 
the  effect  of  the  New  England  seamount  chain  on  the  Gulf  Stream 
(Hurlburt  and  Thompson,  1984). 


2.0  Model  Formulation 
2.1  Model  Equations 


The  vertically  integrated  equations  of  motion  used  in  the  n-layer 
finite  depth  hydrodynamic  model  are,  for  /t=  1.  .  .n: 


^  +  (V  .  Vfc  +  Vfc  .  V)n+kxfVk  = 

Bt 

-/li  I  GkiV{hi-Hi) 

i  =  1 

+  max  (0,o)i)v*  +  1  -  (max  (O.-o)^)  +  max  (O.o)*  _  j))v^ 

+  max  (0,-(o*  _  i)v*  _  1 

+  _  1  -Xi  j/po  +  Anhk  V^n 

r^h  -* 

~  +  V  ■Vk=(jik-(iik-\  0) 

dl 

where 

/i*  =  /t-th  layer  thickness; 

vt  =  ik-th  layer  velocity; 

^k  =hkVk\ 

=  /t-th  layer  thickness  at  rest; 
n  -  1 

//„  =D(x,y)-  Z  Hi; 

i  =  1 

D(n,  >’)  =  total  ocean  depth  at  rest; 
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=  coefficient  of  inierfacal  friction; 

=  coefficient  of  bottom  friction; 
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co^  = 

^  JJ(max(0,a)i)-max(0,(O4)) 

COi 

CO*:  =  ^-th  interface  reference  vertical  mixing  velocity; 

/z*  =  A:-th  layer  thickness  at  which  entrainment  starts;  and 

/i*  =  k-th  layer  thickness  at  which  detrainment  starts. 

This  layered  formulation,  with  transports  V  as  dependent  variables, 
handles  strongly  sloping  topography  especially  well,  as  long  as  it  is 
confined  to  the  lowest  layer.  The  bottom  topography  appears  multipli- 
catively  in  the  pressure  gradient  term,  and  is  differentiated  only  to  the 
extent  that  it  affects  the  velocity  field  in  the  advective  terms.  When 
large-amplitude  topography  is  introduced,  restrictions  on  the  time  step 
and  the  eddy  coefficient  are  affected  only  to  the  extent  that  the  topography 
determines  the  amplitude  of  the  velocity  field.  The  continuity  equation 
is  linear  when  written  in  transport  form,  thus  avoiding  complications 
from  the  nonlinear  advective  term  when  layer  thickness  and  velocity 
are  used  as  variables. 

The  surfacing  of  layer  interfaces  has  always  been  a  problem  for  layer 
models.  The  traditional  solution  has  been  to  make  each  layer  sufficiently 
thick  to  prevent  the  problem  from  occurring.  However,  this  practice 
often  leads  to  a  vertical  structure  that  is  not  in  line  with  that  of  the 
ocean  region  under  study.  The  hydrodynamic  model  is  an  isopycnal 
model  (the  density  is  constant  in  space  and  time  for  each  layer); 
however,  hydromixing  has  been  included  to  prevent  surfacing.  The  physics 
of  hydromixing  are  identical  to  mixing  in  a  layered  thermodynamic 
model  currently  under  development  at  NOARL,  except  that  the  fluid 
entrained  from  another  layer  is  assumed  to  be  at  the  same  density  as  the 
receiving  layer.  It  is  primarily  designed  to  prevent  any  layer  interface 
from  surfacing,  and  depends  only  on  the  layer  thickness  and  three  user- 
supplied  parameters;  a  reference  vertical  mixing  velocity.  0)*;  a  layer 
thickness  at  which  to  start  entrainment,  /z*;  and  a  layer  thickness  at 
which  to  start  detrainment,  /z*. 

In  practice  explicit  detrainment  is  often  deactivated  by  making  /z* 
very  deep.  Mixing  has  been  formulated  to  involve  no  net  transfer  of 
fluid  between  layers,  so  if  there  is  local  entrainment  it  will  be  balanced 
by  region-wide  detrainment.  This  global  mixing  balance  term,  O)*,  has 
been  included  to  allow  long-term  integrations;  however,  it  is  an  additional 
source  of  damping.  For  regions  with  port  forcing  the  layer  thickness 
can  alternately  be  balanced  by  adjusting  the  net  transport  in  or  out  of 
the  region,  thus  removing  the  need  for  the  to*  term.  Hydromixing  is  a 
relatively  new  addition  to  the  ocean  model;  it  significantly  expands  the 
range  of  problems  that  can  be  investigated  with  the  model,  but  it  should 
be  used  only  if  the  .solutions  with  mixing  provide  a  more  realistic 
simulation  than  those  without  mixing. 

A  hydrodynamic  reduced  gravity  model  with  active  layers  has  a 
lowest  layer  that  is  infinitely  deep  and  at  rest,  i.e.,  +  i  =  0,  /z„  +  j  =  «=. 
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and  V/i„  +1  =  0.  The  model  equations  for  the  active  layers  are  identical 
to  those  for  an  «-Iayer  hydrodynamic  finite  depth  model,  except  that 


2.2  Finite  Difference  Equations 


//„  =  constant 

n  —  I  +  t  ~ 

■  1  «(Pn  +  1  -  P;)/Po 

j 

1Q  Po  I  V*  -  V*  +  1 1  {Vi  -  Vi  +  i) 

(0  ^ 

\  max  (0,a)i)-max(0,a)i)-/i*C0i 


ti  = 
and 

(Oi  = 


for  l<  k 
for  /  >  /t 

for  A:  =  0 

for  /:  =  1  .  .  .  n 

for  A:  =  0 

for  A:  =  1  .  .  .  n. 


In  general,  the  reduced  gravity  model  is  more  robust  than  the  finite 
depth  model,  and  an  n  -  1  active-layer  reduced-gravity  model  can  give 
similar  results  to  an  /i-layer  finite  depth  model.  So  it  is  often  cost 
effective  to  develop  a  reduced  gravity  model  of  a  region  first  and  then 
go  on  to  a  finite  depth  model. 


The  model  equations  are  solved  numerically  on  a  staggered  C-grid. 
For  simplicity,  the  P-plane  version  is  described  here.  The  operators  are 
defined  as 


K2(W(z))=imAf 


f  / 
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\  K 


2  +  ■ 


m  A 


-W 


W  (2)mz  =  1/2 


(w 

+  W 

f  mAV 

1  1  2  J 

V  2  )) 

(2) 


where  W  is  a  function  of  the  discrete  variable  r,  A  refers  to  a  space  or 
time  increment,  and  m  is  an  integer  assumed  to  be  1  if  omitted.  Then 
the  finite  difference  form  of  the  finite  depth  hydrodynamic  model  equa¬ 
tions  for  A:  =  1  .  .  .  n  is  given. 

52,(y  =  -d,  (v^uy)+(fvy/j{ 
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for  k  =  0 


//i=  constant.  (14) 


Here,  V  =  (U,V),u  =  U/h^,v  =  V/h^ ,  and  the  layer  index  k  and  the 
centered  time  t  apply  to  all  fields  unless  otherwise  noted. 

For  two  layers  equations  (3)  through  (14)  are  identical  to  those  in 
Hurlburt  and  Thompson  (1980),  except  for  hydromixing  and  the  pres¬ 
sure-gradient  term.  The  original  pressure-gradient  term  involved  the 
free-surfacc  deviation.  t|,  and  the  upper  layer  thickness,  /i],  whereas 
now  the  layer  thickness  deviations,  (/ij  -//j)  and  (hi- Hi),  are  used. 
These  two  forms  are  interchangeable,  since  r\  =  (hi  -  Hi)  -  (hi- Hi) 
and  Hi  is  constant,  but  layer  thickness  deviations  are  prognostic  fields 
in  the  Navy  model;  hence,  it  is  more  natural  to  use  them  directly. 

Only  the  pressure-gradient  term  in  the  momentum  equation  and  the 
mass-divergence  term  in  the  continuity  equation  have  been  made  implicit; 
i.e.,  this  is  a  semi-implicit  formulation  of  the  equations  that  treats  the 
linear  part  of  the  gravity  waves  implicitly.  It  is  unconditionally  stable, 
provided  the  constants  H^  are  such  that  the  explicit  part  of  the  pressure- 
gradient  term  always  acts  to  slow  the  gravity  waves.  The  model  requires 
W*  to  be  greater  than  which  guarantees  unconditional  stability. 

The  linear  CFL  limit  (Courant  et  al.,  1928)  on  the  time  step  is 
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and  the  actual  time  step  is  typically  between  90%  and  95%  of  this  limit. 

Expressing  the  n  continuity  equations,  one  for  each  layer,  as  a  single 
« -element  vector  equation,  we  obtain 
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where  5^  =  6^,6^,,  and  hi  contains  all  the  terms  at  time  levels  t 

and  t-At. 

2 

Now,  a  real  eigenmatrix  X  and  eigenvalues  exist,  such  that 
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2.3  Boundary  Conditions 


The  semi-implicit  continuity  equation  in  real  space  consists  of 
« -coupled  elliptic  partial  differential  equations.  By  converting  to  modal 
space,  these  are  decoupled  into  n  independent  two-dimensional  Helmholtz 
equations.  The  decoupling  gives  rise  to  Helmholtz  equations  because 
Hk  and  G*/  (and,  therefore,  X  and  )  are  constants. 

The  time  differencing  is  leapfrog,  but  an  Heun  predictor-corrector 
difference  (Roache,  1976)  is  used  for  start  and  restart.  Every  so  often, 
typically  every  100  to  300  time  steps,  the  solution  for  at  least  three 
adjacent  time  steps  is  averaged  to  exactly  filter  2Ar  time  splitting,  and 
the  model  is  restarted. 


On  the  C-grid:  zonal  boundary  sections  pass  through  V  grid  points, 
meridional  sections  pass  through  U  grid  points,  and  diagonal  boundaiy 
sections  pass  through  both  U  and  V  points.  In  all  cases  the  boundary  is 
exactly  Ax/2  or  A>72  from  the  nearest  h  grid  point.  So  the  coastline 
is  uniquely  defined  by  a  map  indicating  which  h  grid  nodes  are  sea  and 
which  are  land.  See  the  appendix  for  details  of  the  grid  used  by  the 
model. 

Except  at  inflow  and  outflow  boundary  ports,  the  walls  are  rigid  and 
the  no-slip  condition  is  prescribed  on  the  tangential  flow.  These  bound¬ 
ary  conditions  are  implemented  in  the  model  by  setting  all  velocity  and 
transpon  components  that  lie  on  the  boundary  to  zero,  and  setting  all 
components  that  lie  one-half  grid  spacing  on  the  land  side  of  the  boundary 
equal  to  minus  those  that  are  one-half  grid  spacing  on  the  sea  side  of 
the  boundary.  If  a  (7  land  point  has  U  sea  points  directly  to  its  north  and 
south,  or  if  a  V  land  point  has  V  sea  points  directly  to  its  east  and  west, 
then  the  boundary  must  pass  through  the  land  point.  At  a  port  setup  for 
inflow  only,  the  transport  ,  or  velocity  v^,  are  external  forcing  fields. 
The  flow  at  a  port  setup  for  outflow  is  predicted  in  three  steps. 

First,  a  modified  Orlanski  (1976)  radiation  condition  is  used.  Con¬ 
sider,  for  example,  a  port  on  the  northern  boundary  of  the  region.  Then, 
‘  .  r  each  component  of  transport. 


;  +  A/ 


cAt 

Ay 


I  -  At 


+ 


2cAt  t 


1  + 


cAt 

Ay 


c=  min(A>7A  t.  max(c*,Vi)), 


(19) 


where  B  refers  to  the  outflow  boundary  and  B  -  1  to  one  grid  point 
interior  to  the  boundary,  v*  is  a  reference  outflow  speed  that  is  usually 
set  to  the  average  inflow  speed  at  the  inflow  pons,  and  the  phase  speed, 
Ck  is  determined  locally  for  each  wave  r  dc.  TTiere  are  several  ways 
to  determine  the  local  phase  speed,  c*.  bui  the  simplest  scheme,  with 
c*  =  V* ,  is  often  very  effective.  Overall,  this  boundary  condition  always 
moves  information  from  the  interior  toward,  and  normal  to,  the  open 
boundary. 
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Second,  a  linear  implicit  drag  is  an  option  that  may  be  applied  at  all 
points  on  the  port  where  inflow  is  predicted.  This  option  prevents  the 
development  of  any  unrealistically  large  recirculation  at  the  open  boundary 
by  applying  a  “brake”  on  the  inflow  component  of  the  recirculation. 

The  e-folding  time  scale  of  the  drag  is  chosen  to  be  sufficiently  long 
to  allow  waves  and  eddies  to  migrate  out  of  the  region. 

Third,  the  velocity  across  each  outflow  port  in  a  layer  is  uniformly 
increased  or  decreased  such  that  the  total  inflow  mass  transport,  from 
all  inflow  ports  in  that  layer,  is  matched  by  the  total  outflow  mass 
transport.  This  hard  mass  constraint,  time  step  by  time  step,  on  inflow  and 
outflow  does  not  allow  for  temporary  mass  storage  in  the  region,  and  it 
can  lead  to  problems  when  spinning  up  a  model  from  rest.  However, 
it  allows  the  use  of  excess  mass  as  a  debugging  tool  (since  there  should 
not  be  any  net  change  in  mass  with  this  scheme). 

To  solve  the  Helmholtz  equations  for  /i,  we  require  boundary  values 
for  U  and  V  at  the  new  time  level.  Since  interior  U  and  V  values  are  not 
available  until  after  the  new  h  is  known,  some  boundary  conditions  on 
U  and  V  that  involve  interior  values  at  the  new  time  level  could  not  be 
used  in  a  semi-implicit  model.  However,  simple  boundary  conditions  of 
this  form,  such  as  those  at  a  free -slip  boundary,  do  not  cause  any  dif¬ 
ficulty.  In  the  present  model  at  a  solid,  no-slip  boundary,  U  and  V  are 
both  zero  for  all  time;  at  an  inflow  port,  they  are  prescribed;  and  at  an 
outflow  port,  they  depend  on  interior  values  from  previous  time  steps 
only.  So  in  all  cases  the  required  boundary  values  are  available  when 
the  Helmholtz  equation  is  solved.  Note  that  the  value  of  h  cannot  be 
prescribed  at  an  open  boundary  in  this  model. 

The  allowed  external  forcing  fields  for  the  model  are  bathymetry,  2.4  External  Forcing  Input 
wind  stress,  and  inflow  boundary,  velocity,  or  transport.  All  of  these 
forcing  fields,  except  for  bathymetry,  are  optional  and  can  be  omitted 
from  a  given  model  run. 

Bathymetry  is  defined  on  the  same  grid  as  h,  and  it  is  usually 
read-in  from  a  file  at  the  beginning  of  the  model  run.  It  must  be  sufficiently 
deep  to  be  always  below  the  lowest  layer  interface,  but  is  otherwise 
arbitrary.  However,  it  is  advisable  to  smooth  out  any  2Ax  or  ZAy 
components  that  are  present  before  using  a  bathymetry  field  in  the 
model.  The  bathymetry  field  is  also  used  to  define  the  model  coastline. 

If  a  point  has  a  positive  bathymetric  value,  then  it  is  taken  as  sea; 
otherwise,  it  is  a  land  point.  Any  land-sea  geometry  is  allowed;  there 
is  no  limit  on  the  number  of  islands  in  the  geometry,  and  single-point 
islands  are  allowed.  Since  the  bathymetry  file  also  defines  the  coast¬ 
line.  it  must  be  present.  However,  in  reduced  gravity  models  it  has  no 
other  effect,  and  in  finite-depth  models  the  amplitude  of  the  topography 
can  be  specified.  An  amplitude  of  zero  corresponds  to  a  flat-bottom 
case,  and  an  amplitude  of  one  corresponds  to  the  topography  exactly  as 
input.  Any  amplitude  may  be  used,  but  typically  the  bathymetry  might 
be  defined  to  reach  the  surface,  and  an  amplitude  slightly  less  than  one 
would  be  used  to  ensure  that  the  actual  model  bathymetry  always  remained 
below  the  lowest  interface. 

The  wind  stress  fields  are  defined  on  the  U  and  grids.  Wind  stresses 
are  usually  read-in  from  a  file,  or  files,  during  the  model  run.  The  time 
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3.0  Model  Code 

3.1  Programming  Language 


3.2  Use  of  Third-Party 
Modules 


intervals  between  wind  input  are  arbitrary  and  can  change  from  one 
wind  record  to  the  next.  At  any  given  time  step  the  model  wind  stress  is 
formed  by  a  linear  interpolation  of  the  nearest  input  wind  stress  before 
that  time  and  the  nearest  input  wind  stress  after  that  time.  Since  the 
scale  of  wind  features  is  usually  far  larger  than  the  scale  of  oceanic 
features,  only  one  out  of  every  two  wind  stress  values  is  input  and 
saved  in  each  direction,  i.e.,  only  one-fourth  of  the  points  are  in  the 
wind-stress  file.  See  the  appendix  for  more  details  of  the  wind  grid. 
The  wind-stress  values  at  the  remaining  three-fourths  of  the  points  are 
obtained  by  linear  interpolation  in  space.  Wind  fields  must  be  converted 
to  wind  stresses  at  the  surface  before  they  are  input  to  the  model.  If  a 
more  sophisticated  time  interpolation  scheme  (than  the  linear  interpolation 
used  by  the  model)  is  required,  interpolated  wind  stress  fields  are  included 
in  the  data  file  between  the  times  that  actual  winds  are  available. 

At  an  inflow  boundary  in  any  layer,  the  inflow  angle  and  either  the 
velocity  magnitude  and  profile  or  the  transport  magnitude  and  profile 
can  be  specified.  These  quantities  are  not  allowed  to  change  with  time, 
except  on  restart  when  a  new  inflow  angle  and  new  magnitudes  can  be 
specified,  and  the  model  will  spin  up  from  the  old  to  the  new  values. 
At  the  start  of  the  run  the  magnitudes  spin  up  from  the  initial  state  to 
the  specified  values.  Outflow  boundaries  use  a  radiation  condition  and 
do  not  use  any  specified  data. 


The  ocean  model  is  written  entirely  in  the  American  National  Stan¬ 
dard  Institute’s  standard  FORTRAN  77,  with  the  exception  of  the 
following  extensions: 

•  NAMELIST  I/O  (input/output)  is  used  to  read-in  model  parameters. 

•  DATA  statements  within  normal  subroutines  are  used  to  initialize 
variables  in  some  named  COMMON  areas. 

•  Some  I/O  routines  are  machine-dependent,  and  may  include  non¬ 
standard  FORTRAN  features  on  some  machines,  e.g.,  LOGICAL*!  or 
BYTE  data  types  might  be  used. 

All  the  machine-dependent  code  is  accessed  via  a  single  subroutine. 
The  reason  for  using  machine-dependent  I/O  is  to  allow  the  data  to  be 
machine-independent  and  compact.  The  only  standard  FORTRAN  method 
of  producing  machine-independent  data  is  formatted  I/O,  but  this  procedure 
is  slow  and  wastes  storage  (disk)  space.  With  the  current  scheme,  model 
data  generated  on  one  computer  can  be  used  later  on  another  computer 
of  a  different  type,  either  to  produce  graphical  output  or  to  restart  the 
model  and  continue  to  run. 

LINPACK  and  EISPACK  are  the  only  third-party  modules  called  by 
the  model.  These  modules  solve  simultaneous  linear  equation  and 
eigenvalue  problems,  respectively.  Both  are  in  the  public  domain, 
and  highly  optimized  versions  are  often  available  in  the  standard  math 
library  on  supercomputers. 

The  Helmholtz  equations  generated  by  the  semi-implicit  time  scheme 
are  solved  using  a  locally  written  standard  FORTRAN  77  package  that 
has  been  highly  optimized  for  use  on  vector  supercomputers.  It  uses  a 
fast  direct  solver  for  rectangular  regions  (Swarztrauber,  1977)  and  the 
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capacitance  matrix  technique  for  general  regions  (Buzbee  et  al.,  1971; 

Wallcraft,  1980).  The  fast  direct  solver  includes  an  efficient,  standard 
FORTRAN  77  package  for  fast  Fourier  transforms  of  length  4  x  2^  or 
5  X  2^  or  6  X  2^  (Moore.  1989). 

Graphical  output  is  produced  as  a  postprocessing  step  separate  from 
running  the  ocean  model.  The  National  Center  for  Atmospheric  Research 
graphics  package,  together  with  a  level  2B  implementation  of  the  Graphics 
Kemal  System  (GKS)  is  used  as  the  basis  for  all  graphical  output.  GKS 
is  an  international  standard  for  two-dimensional  graphics  and  is  widely 
available — at  least  six  commercial  implementations  are  available  for 
Digital’s  VAX  computer.  Since  the  model  data  is  machine -independent, 
it  can  be  copied  to  any  suitable  workstation  or  computer  to  produce  the 
graphics. 

At  user-specified  intervals  the  model  checks  that  the  prognostic  fields  3.3  Run  Interruption 
are  within  reasonable  ranges,  and  it  will  terminate  if  they  are  not.  This 
inspection  is  useful  on  such  computers  as  the  Cray,  which  have  very 
large  exponent  ranges  and  take  a  long  time  to  overflow  if  the  model 
blows  up.  Standard  interactive  debuggers  are  not  usually  very  useful 
when  a  model  overflows,  because  the  problem  often  occurs  much  ear¬ 
lier.  The  usual  procedure  in  such  a  case  would  be  to  restart  from  the 
last  available  restan  data  set  and  frequently  write  out  fields  for  off-line 
graphical  display. 

At  another  user-specified  interval  a  complete  data  set  suitable  for 
restarting  the  model  is  output.  An  internal  check  ensures  that  this  interval 
corresponds  to  no  more  than  a  set  amount,  typically  30  minutes,  of 
central  processing  time.  If  the  model  is  interrupted  for  any  reason,  then 
it  can  be  restarted  at  the  time  of  the  latest  restart  data  set.  The  restan 
involves  changing  the  start  date  of  the  run,  specifying  the  correct  files 
to  read  the  restart  data  from,  and  resubmitting  the  job. 

After  a  restart  data  set  is  output,  the  amount  of  central  processing 
time  the  job  has  left  is  checked,  and  the  program  will  terminate  if 
the  time  is  not  sufficient  to  reach  the  next  restart  dump.  This  procedure 
ensures  that  enough  time  is  left  after  the  run  to  save  files,  etc.  This 
process  is  particularly  important  on  computers  that  delete  local  files  at 
the  end  of  a  run. 

The  model  has  been  run  on  the  following  computer  systems:  3.4  Implementation 

•  Cray  X-MP  under  COS,  using  CFT  and  CFT77  compilers. 

•  Cray  2  under  UNICOS,  using  the  CFT77  compiler. 

•  Cray  Y-MP  under  UNICOS.  using  the  CFT77  compiler. 

•  Cyber  205  under  VSOS,  using  the  FTN200  compiler  and  VAST, 

VAST  II,  or  KAP  precompiler.  A  precompiler  is  essential  for  efficiency 
on  the  205. 

•Convex  C-1  under  UNIX,  using  the  fc  compiler. 

•  Alliant  FX/80  under  Concentrix  (UNIX),  using  the  FORTRAN 
compiler. 

•  VAX  (pVAX  II,  1 1/780,  1 1/785,  8650,  and  8800)  under  VMS.  using 
the  VMS  FORTRAN  compiler. 

•  Sun  4/1 10  under  SunOS  (UNIX),  using  Sun  FORTRAN. 

•  Compaq  DESKPRO  386/20  under  Interactive  System’s  386/ix  (UNIX 
System  V/386),  using  Microway’s  NDP  FORTRAN  386. 
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All  the  required  compiler  directives  are  present  in  the  code  for  each 
of  the  different  vectorizing  compilers.  The  model  contains  internal 
documentation  detailing  any  minor  modifications  to  the  code  required 
on  each  machine.  These  modifications  are  usually  concerned  with  OPEN, 
since  it  is  difficult  to  make  this  statement  completely  machine-independent. 
On  almost  all  of  these  machines  the  model  runs  in  single  precision, 
which  contains  32  bits  on  everything  except  the  Grays.  The  semi-implicit 
model  must  run  in  64-bit  mode  on  the  Cyber  205  because  it  has  a 
particularly  inaccurate  implementation  of  32-bit  arithmetic.  The  precision 
of  the  model  is  determined  by  an  IMPLICIT  statement  at  the  head  of 
each  routine. 

Most  of  the  work  required  to  implement  the  model  on  a  new  computer 
system  is  in  writing  the  machine-dependent  routines  for  machine- 
independent  I/O  and  in  porting  the  Helmholtz  equation  package. 
Implementing  the  I/O  routines  on  a  new  machine  requires  knowledge  of 
how  to  efficiently  handle  data  that  are  smaller  than  the  usual  machine 
word  length.  The  existing  Sun  FORTRAN  version  will  work  on  most 
machines  that  run  UNIX,  but  it  would  probably  have  to  be  modified  to 
attain  maximum  efficiency  on  a  new  UNIX-based  vector  or  parallel 
computer.  The  Helmholtz  equation  package  is  written  in  standard 
FORTRAN  77,  and  test  programs  are  available  to  check  that  it  is  producing 
correct  results.  The  package  is  designed  to  be  automatically  vectorized 
by  a  wide  range  of  compilers,  but  additional  optimization  (extra  compiler 
directives,  for  example)  would  probably  be  required  to  attain  full  efficiency 
on  a  new  vector  or  parallel  computer.  The  Helmholtz  equation  package 
is  treated  as  a  “black  box”  by  the  ocean  model,  so  it  could  be  replaced 
by  any  other  Helmholtz  equation  package  that  includes  the  capability  to 
solve  problems  with  homogeneous  Neumann  boundary  conditions 
located  halfway  between  grid  points  on  the  /i-grid. 


4.0  Using  the  Model  The  model  code  contains  internal  documentation  that  fully  describes 

4.1  Overview  the  various  user-specified  model  parameters  and  data  sets.  The  details 

of  this  information  will  not  be  reproduced  here,  since  it  is  likely  to 
change  slightly  with  each  new  release  of  the  model.  The  internal  docu¬ 
mentation  will  always  be  current,  so  this  user’s  guide  should  be  used 
only  to  gain  a  general  idea  of  how  the  model  works.  Some  of  the 
internal  documentation  applies  only  to  the  thermodynamic  version  of 
the  model,  which  is  still  under  development  and  has  not  been  released 
for  general  use. 


4.2  The  Model  Label  Each  experiment  has  a  user-supplied  label  associated  with  it.  This 

label  appears  in  all  data  files  and  on  all  the  standard  graphical  products 
supplied  with  the  model.  The  label  consists  of  a  five-digit  integer  followed 
by  a  four-digit  integer.  The  five-digit  integer,  together  with  the  number 
of  active  layers,  should  uniquely  identify  the  region  and  the  model  used 
for  the  experiment.  It  is  set  by  a  PARAMETER  statement  in  the  model 
code,  so  experiments  run  from  the  same  code  will  definitely  have  the 
same  five-digit  identifier.  The  second,  four-digit  integer  should  uniquely 
identify  the  experiment,  so  it  is  set  at  run  time. 
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The  five  digit  integer  has  the  following  form: 

•  a  two-digit  region  identification  number 

•  a  one-digit  region  version  number 

•  a  one-digit  model  type  number 

-  1  for  explicit  hydrodynamic  models 

-  3  for  semi-implicit  hydrodynamic  models 

•  a  one-digit  model  “geometry”  number 

-  0  for  reduced  gravity  models  on  a  beta-plane 

-  1  for  reduced  gravity  models  on  a  sphere 

-  2  for  finite  depth  models  on  a  beta-plane 

-  3  for  finite  depth  models  on  a  sphere 

The  four-digit  integer  has  the  following  form: 

•  a  three-digit  experiment  identification  number 

•  a  one-digit  subexperiment  identification  number 

The  four-digit  integer  is  divided  by  10  to  give  the  overall  experiment 
number;  for  example,  the  first  experiment  would  probably  be  called 
1.0,  so  the  integer  would  be  0010.  By  convention  each  experiment 
started  from  rest  has  a  subexperiment  identifier  of  0,  and  the  other 
subexperiment  identifiers  are  resers’ed  for  experiments  restarted  from 
that  experiment  with  different  model  parameters  or  different  model  forcing. 

For  example,  at  NOARL  the  region  identification  number  of  11  is 
used  for  world  ocean  regions.  Currently,  six  versions  are  available:  111, 

112,  113,  114,  115,  and  116,  ail  of  which  have  differing  north-south 
extents  and  differing  grid  resolutions.  So  from  an  experiment  label, 
such  as  11611, 0134,  we  can  immediately  tell  that  the  experiment  is  for 
the  world’s  oceans  on  our  latest  grid,  that  the  experiment  is  hydrodynamic, 
and  that  it  uses  spherical  coordinates.  Since  the  experiment  is  13.4,  we 
also  know  that  it  is  a  variation  of  experiment  13.0;  for  example,  13.0 
might  have  been  run  for  50  years  from  rest  using  Hcllermann-Rosenstein 
(1983)  climatological  winds,  and  experiment  13.4  might  have  been 
restarted  from  the  50th  year  of  13.0  with  some  other  wind  set. 

The  following  steps  are  taken  to  set  up  the  model  on  a  new  ocean  4.3  Model  Setup 
region.  4.3.1  On  a  New  Ocean  Basin 

•  Select  the  model  grid  spacing,  the  region  extent,  and  the  initial 
vertical  structure.  These  areas  define  what  kind  of  model  to  use  (spherical 
versus  beta-plane,  reduced  gravity  versus  finite  depth,  grid  si^c,  number 
of  layers). 

•  Always  start  from  the  Gulf  of  Mexico  code  as  originally  distributed 
(with  LABEL1=990). 

•  Modify  the  PARAMETERS  IH,  JH,  AND  KH  throughout. 

•  Define  the  PARAMETER  LABEL  1  to  uniquely  identify  the  region/ 
model  combination  (see  above). 

•  Create  a  new  topography  file  and  possibly  modify  PARAMETERS 
LXX,  IXX,  lYY  in  the  subroutine  XXLND.  ' 

•  Modify  many  PARAMETERS  in  subroutine  CHSLVX  or  SHSLVX 
to  set  up  the  Helmholtz  equation  .solver  for  the  new  region,  and  modify 
all  the  entries  in  MHSLV  to  call  cither  CHSLV  or  SHSLV  as  required. 

•Modify  XXPOR  (ZPlOl,  ZPOOl,  etc.)  to  implement  the  required 
ports.  Use  a  dummy  version  of  XXPOR  if  there  is  no  port  forcing. 

•  Possibly  define  a  new  analytic  wind  stress  function  in  subroutine 
XXAWND. 

•  Possibly  generate  one  or  more  wind  stress  files. 


« 
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4.3.2  On  a  New  Type 
of  Computer 


5.0  Summary 


6.0  References 
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The  following  steps  are  taken  to  implement  the  model  code  on  a  new 
computer; 

•  Use  a  FORTRAN  77  compiler. 

•  Use  the  NAMELIST  statement  to  read-in  model  parameters;  if 
NAMELIST  is  not  available,  use  READ  (s,*). 

•  Develop  versions  of  the  history  file  I/O  routines  (XHHIS,  ZHHIS, 
ZHIOXX,  ZHTXTOX)  that  read  and  write  history  fields  in  exactly  the 
correct  manner,  i.e.,  as  24-bit,  2’s  complement  integers. 

•  Reduce  the  PARAMETER  NHTMAX  in  subroutine  XHHIS  appro¬ 
priately  if  the  history  file  unit  numbers  must  be  limited  to  less  than 
65  .  .  .  80. 

•  Test  against  an  existing  benchmark  run  on  another  machine  if  pos¬ 
sible.  Otherwise,  check  that  mass  is  conserved  by  running  in  single  and 
double  precision  and  comparing  “mean  height  anomalies,”  which  should 
be  smaller  by  a  factor  of  about  10*  in  double  precision.  The  U,  V’,  and 
h  fields  should  closely  agree  in  the  two  cases. 


This  report  is  a  users  guide  to  the  Navy’s  hydrodynamic  (isopycnal) 
nonlinear  primitive  equation  layered  ocean  circulation  model.  Tlie  model 
retains  the  free  surface  and  uses  a  semi-implicit  time  scheme  that  treats 
all  gravity  waves  implicitly.  It  can  handle  full-scale  bottom  topography, 
provided  it  is  confined  to  the  lowest  layer,  and  an  arbitrary  coastline 
geometry. 

The  model  has  been  in  use  at  NOARL  for  more  than  10  years  for 
simulations  of  the  ocean  circulation  in  the  Gulf  of  Mexico,  The  Caribbean 
Sea,  the  Alboran  Sea,  the  western  Mediterranean  Sea,  the  Indian  Ocean, 
the  western  North  Atlantic,  the  North  Pacific  Ocean,  and  the  global 
ocean.  The  model  code  is  being  made  available  to  the  ocean  modeling 
community  in  conjunction  with  the  issuing  of  this  report. 

The  vertically  integrated  equations  of  motion  and  their  finite  differ¬ 
ence  discretization  on  a  C-grid  have  been  presented,  as  has  a  description 
of  the  semi-implicit  time  scheme,  the  boundary  conditions,  and  the 
external  forcing. 

The  model  code  contains  internal  documentation  that  fully  describes 
the  user  specified  model  parameters  and  data  sets.  This  report  also 
contains  general  information  about  how  to  use  the  model,  in  particular, 
how  to  set  it  up  for  a  new  ocean  region  and  how  to  port  it  to  a  new 
computer  system. 
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Appendix 
The  C-Grid 


The  following  figures  show  the  relative  position  of  the  U,  V  and  h 
staggered  grids  in  the  C-grid  scheme  used  by  the  ocean  model.  They 
also  show  how  grid  points  are  numbered  within  two-dimensional  arrays 
in  the  model.  The  size  of  the  mesh  array  is  chosen  to  just  cover  the 
model  region,  but  because  the  array  is  rectangular  a  large  fraction  of 
the  array  might  represent  “land”  for  some  regions.  In  this  case,  the 
region  is  a  simple  rectangle  12Ax  by  llAy,  so  all  mesh  arrays  are 
dimensioned  13  by  12  (i.e.,  IH  =  13  and  JH  =  12),  and  very  little  of  the 
array  is  “land.”  The  sides  of  the  region  are  on  the  U  grid,  at  I  =  1  and 
1=12,  and  the  top  and  bottom  are  on  the  V  grid  at  J  =  1  and  J  =  1 1.  If 
the  region  has  solid  walls,  these  U  and  V  points  will  be  set  to  zero  as 
part  of  the  no-slip  boundary  condition.  On  the  h  grid  all  the  nodes  in 
rows  J  =  1  and  J  =  12  and  in  columns  1  =  1  and  1=13  are  outside  the 
region.  They  are  never  used  by  the  model,  but  would  be  set  to  negative 
values  in  the  topography  array  to  indicate  that  they  are  over  land.  All 
other  nodes  in  the  topography  array  would  have  positive  values 
representing  the  ocean  depth  at  that  point.  On  the  U  grid  all  nodes  in 
column  1  =  13  are  outside  the  region  and  are  never  used.  All  nodes 
in  rows  J  =  1  and  J  =  12  are  one-half  grid  point  outside  the  region  and, 
assuming  the  boundary  is  a  solid  wall,  these  will  be  set  equal  to  minus 
the  corresponding  nodes  in  rows  J  =  2  and  J  =  1 1  as  part  of  the  no-slip 
boundary  condition.  On  the  V  grid  all  nodes  in  row  J  =  1 2  are  outside 
the  region  and  are  never  used.  All  nodes  in  columns  I  =  1  and  1  =  13  are 
one-half  grid  point  outside  the  region,  and  assuming  the  boundary  is  a 
solid  wall,  these  will  be  set  equal  to  minus  the  corresponding  nodes  in 
columns  I  =  2  and  I  =  12  as  part  of  the  no-slip  boundary  condition. 

Wind  stress  components  are  stored  only  at  every  second  point,  as  is 
indicated  in  the  second  figure.  Wind  stresses  must  be  defined  at  all 
points  on  these  subsampled  grids  because  the  values  over  land  may  be 
used  in  the  interpolation  for  a  sea  point. 
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IH  =  13  and  JH  =  12 
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IH  =  13  and  JH  =  12,  IW  =  7  and  JW  =  6 


Staggered  wind  grid. 
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The  wind  gridpoints  are  u,  v  and  the  nonwind 
gridpoints  are  x,  y. 
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